A little context

This pandemic context, in which a completely new virus is spreading among the world’s population, presents many unknowns and uncertainties.

From the very beginning of the epidemic, one of the main concern about this emerging disease is its fatality rate.

This information is important as it allows to quickly get a sense of the most dramatic consequences that this disease can potentially have in the event of a widespread pandemic.

Very soon after the outbreak began, we began to receive data, notably from China, about the number of patients identified (tested) and the number of victims of the new virus SARS-CoV-2. The current consensal source for the collection, organisation and synthesis of this data can be found here: Johns Hopkins University Dashboard

The legendary formula

Seen like that, this formula seems to be a bit vague (also note that it can be written in different forms). We won’t decipher it in detail here, we will simply note some elements needed to shed some light on why this formula is useful:

• P(A) means probability of A
• |  means knowing
• P(A|B) means probability of A knowing B

From these few details, we can see that the formula gives us the overall probability of A knowing B on the basis of the probability of A, the probability of B and the probability of B knowing A. This is the reason why Laplace spoke of inverse probabilities ;)

If we replace A by fatality rate and B by data we see that our formula allows us to know the probability of the fatality rate knowing our data on the basis of the probability of the data knowing the fatality rate.

This may give us the impression of a vicious circle but in reality it is not. Let’s look at the following illustration, which is a readapted version of the formula:

This version of the formula simply tells us that the probability of a fatality rate knowing our data is proportional to the product of the probability of the data knowing that rate and the probability of the rate.

This formula therefore simply indicates that combining our prior on the fatality rate with the probability that our data can be generated in a world where such a rate is true allows us to calculate the probability of this rate knowing our data, i.e. our posterior.

Implementing a generative model

Of course we have a very general description of the principle of the generative model, but in practice how do we go about it?

Let’s take an example, illustrate it and implement it to get a better insight.

Let’s take the scenario of a MERS-CoV epidemic outbreak. Now imagine that we want to simulate a situation where an epidemiologist goes to this hot spot and identifies 160 infected people. What number of deaths could we expect? Of course we could simply calculate 30% of 160 patients (i.e. 48 victims) and then be satisfied with the result.

But we know that 30% is a rough estimate of the proportion of deaths around which our simulated data should fall, but, because of the randomness, it is not an exact and absolute value at all times. We will therefore use a probability distribution whose parameter will be 30% and which will allow us to produce a number of victims for 160 virtual patients.

Our generative model will therefore implement this probability distribution. In our case, the output of the generative model must tell us for each of the 160 virtual patients whether they are victims or not. The probability distribution capable of meeting these criteria is the binomial distribution.

The following diagram illustrates the process we will implement to generate a simulated data set

In practice, the R code for this particular case of a generative model is as follows:

victimes_virtuelles <- rbinom(1, size = 160, prob = 0.30)
victimes_virtuelles

## [1] 52

If you execute this line of code several times in a row you will find that with each new execution the number of virtual victims changes and orbits more or less around the value 48.

We have therefore been able to simulate a number of virtual victims based on a given fatality rate.

Go backwards to the parameter value

Our problem is that what we’re interested in is the exact opposite ! We have data and we are looking to go backwards to the case-fatality rate parameter, more technically speaking, we are looking to make a Bayesian inference.

This little video from the Mental Hygiene channel exposes this same problem with dice throws. We know that a dice was rolled and it produced a 6. What we are trying to find out is: what type of dice was rolled? Is it a 4, 6, 8, 12, or even 20 sided dice? As with our case-fatality rate, Bayesian inference is required to determine which dice was most likely rolled.

We will find that despite this apparent difficulty, we are on the right track to move from Prior to Posterior.

Implementation of the Prior

In practice, the 10,000,000 fatality rate values are first drawn down

We have defined an uninformed Prior, for which each case-fatality rate can be equiprobably between 0 and 100%, we can implement it as follows:

nb_tirages <- 10000000
prior <- runif(n=nb_tirages, min=0, max=1)

In this line of code, we create a list called prior in which we place the results of the 10000000 prints of values between 0 and 1 (1 equals 100%).

We can illustrate it graphically as follows:

hist(prior, xlim=c(0,1), ylim=c(0,600000), 
     main="A Priori", xlab="Taux de létalité", 
     ylab="Fréquence", col="grey")

We can see from this histogram that each case-fatality rate value is drawn in almost equal quantities.

Generative Model implementation

Then, all these rate values are supplied to the generative model.

Let’s first implement this one

modele_generatif <- function(taux, nb_cas) {
  nb_morts <- rbinom(1, size = nb_cas, prob = taux)
  nb_morts
}

This very simple function takes as input the fatality rate and the number of cases for which death data is to be generated.

As an output, it returns the number of simulated deaths.

Data simulation

We can now provide all our fatality rates to the generative model to simulate the data:

data_cas = 271
morts_simu <- rep(NA, nb_tirages)
for(i in 1:nb_tirages) {
  morts_simu[i] <- modele_generatif(taux = prior[i], nb_cas = data_cas)
}

This small piece of code produces the list of the 1000000 simulated data.

We will then be able to filter this data in order to only keep the ones that correspond to our real data (i.e. 1 death).

posterior <- prior[morts_simu == 1]

We obtain a list of all fatality rate values that among the 10,000,000 equiprobable values from our Prior distribution, when supplied to the generative model, resulted in data that matched our fatality data for Day 1.

We must check that the number of simulated data matching our actual data is sufficient.

length(posterior)

## [1] 36550

This list is none other than our posteror rate distribution. We can represent it by the following histogram

hist(posterior,  breaks=40, axes=F,
     main="A posteriori", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))
axis(2)
axis(1, at=seq(0, 0.05, 0.001))

From this distribution, we can make an estimate of the case-fatality rate interval in which 95% of the drawn case-fatality rates that resulted in simulated data compatible with our real data.

This is the 95% Credible Interval

quantile(posterior, c(0.025, 0.5, 0.975))

##         2.5%          50%        97.5% 
## 0.0008818943 0.0061149305 0.0198645343

mean(posterior)

## [1] 0.00725743

It can hence be seen that 95% of the rates that produced compatible data are in the range of 0.08% and 2% with a median rate of 0.6% and an average rate of 0.7%.

The impact of our data has therefore visibly affected our potential vision of the fatality rate of the disease. We have moved from a situation where this rate could equiprobably take any value between 0 and 100% to a possible value between 0.08% and 2% considering our data at day 1.

The following chart illustrates our transition from the prior to the posterior (note that the scales of the two histograms are very different here).

Prudence

This prediction is interesting, but isn’t it a bit risky to make predictions when our estimates are based only on taking into account information from 271 identified cases and 1 death?

Of course the answer is yes, but let us point out that the approach is still interesting and that it allows us to produce anticipations with little data and, above all, a certain quantification of uncertainty.

Before starting to take into account the data of the following days, let’s dwell for a while on the effect of the data we take into account on the Posterior.

Village A

Let us take again the code executed previously by replacing the data of the first day in Germany by the data of the village A. Here all the stages of simulation are the same (morts_simu) only the filtering of the simulated data is different: where we kept only the simulated data for which the number of deaths was 1 now we will keep only the simulated data for which the number of deaths is 135.

posterior_vill_A <- prior[morts_simu == 135]

length(posterior_vill_A)

## [1] 36505

hist(posterior_vill_A,  breaks=40, axes=T, xlim=c(0, 1), 
     main="A posteriori Village A", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.4))

quantile(posterior_vill_A, c(0.025, 0.5, 0.975))

##      2.5%       50%     97.5% 
## 0.4382669 0.4981636 0.5574826

mean(posterior_vill_A)

## [1] 0.4981637

As you can see, this time the rate estimated in the Posterior ranges between 43% and 55%. We therefore confirm that for a draw of the same number of patients, the number of deaths correctly influences the estimation of the Posterior rate.

As a reminder, here is the histogram of the distribution of the case-fatality rate for first-day data in Germany

hist(posterior,  breaks=40, axes=T, xlim=c(0, 1), 
     main="A posteriori (prior uniforme)", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))

Let’s look at the data for village B

Village B

This time the data at our disposal is much less. Let’s look at the code again to see the effect of this particular situation.

We need to generate new simulated data this time by drawing 6 patients for each rate of the prior distribution.

data_cas = 6
morts_simu_B <- rep(NA, nb_tirages)
for(i in 1:nb_tirages) {
  morts_simu_B[i] <- modele_generatif(taux = prior[i], nb_cas = data_cas)
}

We filter the compatible simulated data (3 deaths)

posterior_vill_B <- prior[morts_simu_B == 3]

length(posterior_vill_B)

## [1] 1425887

Let’s look at the histogram of the Posterior rate distribution.

hist(posterior_vill_B,  breaks=50, axes=T, xlim=c(0, 1), 
     main="A posteriori Village B", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))

Compared to the village A our distribution of probable rates is still centred around 50% however, the spread of the distribution is much wider.

quantile(posterior_vill_B, c(0.025, 0.5, 0.975))

##      2.5%       50%     97.5% 
## 0.1843207 0.5002218 0.8164351

mean(posterior_vill_B)

## [1] 0.5002084

This time the Posterior case-fatality rate ranges possibly between 18% and 81%. We interpret this increase in the 95% credence interval as a measure of our increased uncertainty about the actual case-fatality rate.

As a final comparison, let’s examine what we obtain with the data from village C

Village C

data_cas = 10000
morts_simu_C <- rep(NA, nb_tirages)
for(i in 1:nb_tirages) {
  morts_simu_C[i] <- modele_generatif(taux = prior[i], nb_cas = data_cas)
}

We filter out compatible simulated data (5000 deaths)

posterior_vill_C <- prior[morts_simu_C == 5000]

length(posterior_vill_C)

## [1] 1035

Histogram of the distribution of posterior rates.

hist(posterior_vill_C,  breaks=25, axes=T, xlim=c(0, 1),
     main="A posteriori Village C", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0,1,0,0.5))

quantile(posterior_vill_C, c(0.025, 0.5, 0.975))

##      2.5%       50%     97.5% 
## 0.4908049 0.4998180 0.5102305

mean(posterior_vill_C)

## [1] 0.5000329

We observe here that very large amounts of data have the effect of reducing the credence interval of the most probable rates. In other words, the more data we incorporate into our analysis, the less uncertainty we will have.

This is one of the great strengths of Bayesian data analysis, not only does it allow estimates to be made from even small amounts of variable data, but it also provides information on the uncertainty associated with that estimate.

Turning Posterior into Prior

When we defined our first Prior we decided that it would range between 0 and 100% with an equiprobability for each rate. Mathematically it is a so-called uniform probability distribution.

This distribution is purely mathematical, in fact when implementing the Prior in the ABC algorithm, we generated a list of 10,000,000 rate values respecting this mathematical law. Our real Prior was thus a mathematical law, which we have materialized through the creation of this huge sample.

At the end of the execution of our algorithm we obtained a posterior made up of the list of rates that had allowed the generative model to simulate identical data to our real data. What we obtained and illustrated through histograms is a sample of rate values. Although for teaching reasons we have qualified this sample as posterior. In reality it is not the real posterior. Indeed, the true posterior is, like the prior, a law of probability. A mathematical object that we have approximated through our Posterior sample.

If this posterior sample is large enough, we can consider that it reflects the law of probability that we were looking for.

In order to re-use our Posterior sample as a Prior we will have to find out, from this posterior sample, what is it’s underlying probability. We can already see that it is no longer a uniform law. The ideal would be to have a relatively flexible law and there happens to be one which might be suitable: the Beta distribution.

The beta distribution has 2 parameters α and β whose values allow to adapt the shape of the distribution.

Another advantage of this distribution is that it is bounded between 0 and 1.

Let’s look at how to create a sample of 1000000 values according to this probability law and make the parameters α and β vary to see its flexibility.

α = 1 and β = 1

beta_sample <- rbeta(n = 1000000, shape1 = 1, shape2 = 1)

hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 1 et β = 1", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

When α = 1 and β = 1 the Beta distribution takes the form of a uniform distribution

Let’s try some other values

• α = 1 and β = 5

beta_sample <- rbeta(n = 1000000, shape1 = 1, shape2 = 5)
hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 1 et β = 5", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

• α = 45 and β = 55

beta_sample <- rbeta(n = 1000000, shape1 = 45, shape2 = 55)
hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 45 et β = 55", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

• α = 7 and β = 3

beta_sample <- rbeta(n = 1000000, shape1 = 7, shape2 = 3)
hist(beta_sample, xlim=c(0, 1),
     main="Prior α = 7 et β = 3", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

• α = 50 and β = 50

beta_sample <- rbeta(n = 1000000, shape1 = 50, shape2 = 50)
hist(beta_sample, xlim=c(0, 1), 
     main="Prior α = 50 et β = 50", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(0.7,0.7,0.7,0.5))

We thus have at our disposal a law of probability whose coefficients can capture the distribution of the posterior rate distribution that we have obtained. However, we now need to retrieve the coefficient values α and β so that we can take our Posterior and set it as an updated Prior.

To do this we will use an estimation function to approximate the parameters of the beta law to our posterior distribution and retrieve their values.

This function is provided by the EnvStats package. The syntax is as follows:

library(EnvStats)

## 
## Attaching package: 'EnvStats'

## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm

## The following object is masked from 'package:base':
## 
##     print.default

beta_model <- ebeta(beta_sample, method = "mle")
coeff_alpha <- beta_model$parameters[[1]]
coeff_beta <- beta_model$parameters[[2]]

In this example we try to find the coefficients of the last beta distribution we tested: α = 50 and β = 50

We obtain the following coefficients

round(coeff_alpha)

## [1] 50

round(coeff_beta)

## [1] 50

Let’s see which coefficients of the beta distribution correspond to our day 1 posterior fatality-rate sample.

post_beta_model <- ebeta(posterior, method = "mle")
coeff_alpha <- round(post_beta_model$parameters[[1]])
coeff_beta <- round(post_beta_model$parameters[[2]])
coeff_alpha

## [1] 2

coeff_beta

## [1] 278

Let’s check whether the probability distribution matches the one we obtained in our posterior sample.

new_prior <- rbeta(n = 10000000, shape1 = coeff_alpha, shape2 = coeff_beta)
hist(new_prior, breaks = 60, axes=F, 
     main="Prior mis à jour", xlab="Taux de létalité", 
     ylab="Fréquence", col=rgb(1,0,0,0.5))
axis(2)
axis(1, at=seq(0, 0.05, 0.001))

In performing this check, we created a new sample based on the beta-distribution for our new Prior.

ABC’s not enough, let’s move on to JAGS…

We could start again as we have done so far, however, if you have tried to implement the previous lines of code, you may have noticed that it takes a long time to execute, even too long if your computer lacks of resources.

The ABC algorithm is computationally inefficient, from now on we will use more effective and therefore faster methods.

We will not dwell here on how these algorithms work, just be aware that although they work differently, they follow exactly the same principles as what we have done so far.

R is not always efficient for performing the simulations we need, however, it is very convenient for manipulating and visualizing our data. We will therefore use another more suitable language that we will call directly through R.

There are many languages suitable for Bayesian data analysis. They are commonly called probabilistic programming languages. Among the best known are: STAN, JAGS, WEBPPL, PYMC and many others.

We will be using JAGS but be advised that any other language will do just as well.

install.packages(rjags)

First update

First we need to implement the model by specifying the Prior and the law that rules the generative model. We have established our Prior distribution in the previous sections and we know that the distribution of our generative model is a binomial distribution.

The implementation consists of three stages:

• The definition of the model
• The compilation of the model
• The simulation of the data

library(rjags)

#DEFINE the model
covid_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

# Real Data
cases <- 802
deaths <- 2

#COMPILE the model
covid_jags <- jags.model(textConnection(covid_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
covid_sim <- coda.samples(model=covid_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

Visualization of the posterior

We can therefore go straight to the point and visualize the distribution of Posterior rates resulting from our analysis by integrating the data from day 2.

For a better reading we have added the probability distribution of the prior in red and then superimposed the distribution of the posterior in green.

# Create a df with covid_sim chains and another with the prior
covid_chains<- data.frame(covid_sim[[1]], iter = 1:10000)
df_prior_2 <- data.frame(new_prior)
colnames(df_prior_2)=c("p")

#PLOT the simulated posterior
ggplot(covid_chains, aes(x = p)) + 
  geom_density(data=df_prior_2, fill="red", alpha=0.4) +
  geom_density(fill="green", alpha=0.4) +
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Prior vs posterior jour 2')+
  theme(panel.background = element_blank())

As you can see, these new data have changed the probability distribution of the case-fatality rate.

Let’s see what the new 95% credence interval is.

quantile(covid_chains$p, c(0.025, 0.5, 0.975))

##         2.5%          50%        97.5% 
## 0.0009761953 0.0033984733 0.0081180253

We therefore now believe that the case-fatality rate per case is probably between 0.09% and 0.8%.

We have indeed carried out a Bayesian update of our a priori.

The effect of the Prior

We have previously spent some time examining the effect of the data on the posterior in Bayesian inference. But we have not yet taken the time to study the Prior’s effect. Now that we have updated our Posterior based on a new Prior, we need to better understand what the effect of the Prior is.

We will thus resume our fictive example with a village D in which the data are 50 cases and 25 deaths (50% case-fatality rate). This time, instead of varying the data, we will vary the Prior.

For each graph we display in red the Prior used, in green the obtained Posterior and in grey the previous Posteriors obtained with the previous Priors tested.

With a non informative prior

• α = 1 and β = 1

coeff_alpha <- 1
coeff_beta <- 1
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)

#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(fill="green", alpha=0.4)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior non informatif')+
  theme(panel.background = element_blank())

old_post_vill_D <- vill_D_chains

Let’s repeat this analysis by taking different priors:

• α = 1 and β = 5

coeff_alpha <- 1
coeff_beta <- 5
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")


#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(data=old_post_vill_D, fill="grey", alpha=0.4) +
  geom_density(fill="green", alpha=0.4)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior α = 1 et β = 5')+
  theme(panel.background = element_blank())

old_post_vill_D_2 <- vill_D_chains

• α = 45 and β = 55

coeff_alpha <- 45
coeff_beta <- 55
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")


#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(data=old_post_vill_D, fill="grey", alpha=0.2) +
  geom_density(data=old_post_vill_D_2, fill="grey", alpha=0.4) +
  geom_density(fill="green", alpha=0.4)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior α = 45 et β = 55')+
  theme(panel.background = element_blank())

old_post_vill_D_3 <- vill_D_chains

• α = 7 and β = 3

coeff_alpha <- 7
coeff_beta <- 3
cases <- 50
deaths <- 25


Vill_D_prior_sample <- rbeta(n = 1000000, shape1 = coeff_alpha, shape2 = coeff_beta)
vill_D_prior <- data.frame(Vill_D_prior_sample)
colnames(vill_D_prior)=c("p")


#DEFINE the model
vill_D_model <- "model{
         # Likelihood model for X
         X ~ dbin(p, n)

         # Prior model for p
         p ~ dbeta(a, b)
        }"

#COMPILE the model
vill_D_jags <- jags.model(textConnection(vill_D_model), data=list(a=coeff_alpha, b=coeff_beta, X=deaths, n=cases),
                          inits=list(.RNG.name="base::Wichmann-Hill", .RNG.seed=100))

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 1
##    Unobserved stochastic nodes: 1
##    Total graph size: 5
## 
## Initializing model

#SIMULATE the posterior
vill_D_sim <- coda.samples(model=vill_D_jags, 
                         variable.names=c("p"), 
                         n.iter= 10000)

# Create a df with covid_sim chains and another with the prior
vill_D_chains<- data.frame(vill_D_sim[[1]], iter = 1:10000)

#PLOT the simulated posterior
ggplot(vill_D_chains, aes(x = p)) + 
  geom_density(data=vill_D_prior, fill="red", alpha=0.4) +
  geom_density(data=old_post_vill_D, fill="grey", alpha=0.2) +
  geom_density(data=old_post_vill_D_2, fill="grey", alpha=0.4) +
  geom_density(data=old_post_vill_D_3, fill="grey", alpha=0.6) +
  geom_density(fill="green", alpha=0.8)+
  xlab('Taux de létalité') +
  ylab('Densité de probabilité') +
  ggtitle('Posterior avec un prior α = 7 et β = 3')+
  theme(panel.background = element_blank())

We see here that the same data do not lead to the same Posterior depending on the Prior we use. Indeed, the Posterior we obtain is a form of compromise between the fatality rate of the data and the Prior distribution.

As for the variation in the amount of data, the more certain the prior is (not widely distributed), the stronger influence of the Prior will be. That is to say that the posterior obtained will be closer to it. Thus in a situation where the Prior has a very low uncertainty and the data are present in small quantities, the Posterior will be very close to the Prior. On the other hand, if the Prior is very uncertain (e.g. non informative prior) and the data are present in large quantities, then the Posterior will be closer to the central tendency of the data.

Let’s now pursue our initial analysis.

Iterative algorithm

To do this we will set up a loop that automatically repeats all the steps we have seen so far for each row of the data table.

To better visualize the result of our algorithm, we will produce for each day a graph of the Posterior distribution.

First, let’s start from the beginning by setting the alpha and beta coefficients of the a priori distribution to the value 1 so that we start from a non informative priori. We take the opportunity to indicate on which data table we want to iterate.

We also build a data table in which we will store the credence interval and median of each posterior distribution.

## Empty df for fatality rate update tracing
rate_df = data.frame(day=integer(), Q25=double(), Q50=double(), Q97=double())

## Bayesian model preparation
# select the data table used to feed the process
data_table <- german_data

# Set initial Beta distribution coeffs correponding to a uniform distribution
coeff_alpha = 1
coeff_beta = 1

Then, let’s implement the loop in charge of the iterative calculation of the model.

Here is a descriptive list of the algorithm’s steps

1. Retrieve for each row of the browsed data table, the number of cases identified and the number of deaths recorded.
2. Define the model for JAGS
3. Compile the model
4. Simulate 1000000 parameters of the posterior (Markov Chain Monte Carlo)
5. Store simulations chain in a data frame
6. Estimate the new α and β parameters of the beta distribution for the next iteration.
7. Calculate and add the credence interval as well as the median of the posterior to a table.
8. Create a graph showing the new distribution of the posterior

At the end of the loop we add a last step to create the graph of the evolution of the posterior over the days.

## Bayesian updating 
# For loop that iterates on the german-data table and proceed to the modelling steps
# update the beta distribution to set the daily new priors 
# and plot everyday death rate parameter probability distribution
for (row in 1:nrow(data_table)) {
  
  # 1. retrieve the cases and deaths daily quantities
  cases <- data_table[row, "cases"]
  deaths  <- data_table[row, "deaths"]
  
  # 2. Define the model
  covid_model <- "model{
              # Likelihood model for X
              X ~ dbin(p, n)

              # Prior model for p
              p ~ dbeta(a, b)
              }"
  
  # 3. Compile the model
  covid_jags <- jags.model(textConnection(covid_model), 
                           data=list(a=coeff_alpha, 
                                     b=coeff_beta, 
                                     X=deaths, 
                                     n=cases),
                           inits=list(.RNG.name="base::Wichmann-Hill",
                                      .RNG.seed=100))
  
  # 4. Simulate the posterior
  covid_sim <- coda.samples(model=covid_jags, 
                            variable.names=c("p"), 
                            n.iter= 1000000)
  
  # 5. Create a df with covid_sim chain
  covid_chain<- data.frame(covid_sim[[1]], iter = 1:1000000)

  
  # 6. Approximate the beta parameters  
  # (EnvStats package function)
  coeffs <-  ebeta(covid_chain$p, method = "mle")
  coeff_alpha <- coeffs$parameters[[1]]
  coeff_beta <- coeffs$parameters[[2]]
  

  # 7. Append posterior distribution 2.5%, 50% & 97.5% quantile values to
  # to fatality rate estimation table
  day_number=row
  chain_quant <- quantile(covid_chain$p, c(0.025, 0.5, 0.975))
  Q2 <- chain_quant[[1]]
  Q5 <- chain_quant[[2]]
  Q9 <- chain_quant[[3]]
  quantile_row <- cbind(day_number, Q2, Q5, Q9)
  rate_df <- rbind(rate_df, quantile_row)
  
  # 8. PLOT the simulated posterior
  day <- paste("Jour ", row)
  cas <- paste(cases, "malades")
  morts <- paste(deaths, "morts")
  a <- paste("alpha : ", coeff_alpha)
  b <- paste("beta : ", coeff_beta)
  plot <- ggplot(covid_chain, aes(x = p)) + 
              geom_density(fill="green", alpha=0.4, position = "stack")+
              xlim(0, 0.05) +
              ylim(0, 6000) +
              xlab('Taux de létalité') +
              ylab('Densité') +
              ggtitle(paste0(day, ", ",cas, ", ",morts, ", ",a, ", ",b))
  print(plot)
  file_name <- paste0(path, day, ".png")
  ggsave(filename = file_name)
  }
# End of the loop


# 9. Plot the posterior evolution over daily iterations
ggplot(rate_df, aes(x=rate_df$day, y=rate_df$Q5))+
  geom_line(aes(x=rate_df$day, y=rate_df$Q9), color="grey")+
  geom_line(aes(x=rate_df$day, y=rate_df$Q2), color="grey")+
  geom_ribbon(data=rate_df, 
              aes(ymin=rate_df$Q2,ymax=rate_df$Q9), 
                  fill="lightgrey", alpha="0.7")+
  geom_line()+
  ylim(0, 0.036) +
  xlab('Jours') +
  ylab('Taux de létalité') +
  ggtitle('Évolution de la distribution posterior du taux de létalité')+
  theme(panel.background = element_blank())
ggsave(filename = paste0(path, 'rate_evolution.png'))

Output Graphs

Here is an animation of some of the graphics produced by the execution of the algorithm.

And the graph showing the evolution of the posterior estimate of the case-fatality rate.

Interpretation

We can make two main observations on the basis of these two graphical representations:

• In the first few days, the uncertainty regarding the case-fatality rate is relatively high but decreases over time.

• In the following days, the uncertainty remains globally of a similar amplitude but the rate varies in the direction of a regular increase.

The first finding is largely explained by the fact that the number of cases and the number of deaths remain low in the first few days.

The second is more problematic because it seems that our procedure fails to converge towards a given fatality rate. One of the explanations could be related to the fact that there is a delay between the diagnosis of the disease and the death of the patients (an element that we voluntarily left aside at the beginning of our analysis). This delay implies that the figures obtained at the beginning of the growth phase of the epidemic contain a higher proportion of new patients than of new deaths. Since all patients are new ones, the time for a fraction of them to die has not yet elapsed. The opposite occurs during the declining phase of the epidemic. As the epidemic ebbs and flows backwards, the number of new patients begins to decrease but is not immediately followed by the number of deaths, which may temporarily continue to grow or stabilize. At the end of the epidemic wave there are no new cases, but a continuation of deaths.

Our bias to ignore this time frame has therefore been detrimental to us here. However, we can easily imagine situations in which such a delay does not exist and in which our update strategy would be valid, such situations where individuals would be led to make a yes / no decision for example.